Enseigner avec R

AFROMPA’R Bouaké 2025

INTRODUCTION

On se propoe de passer en revue une partie du programme de l’été à travers un exercice pédagogique associant statistique, cartographie, analyse spatiale, modélisation … Cet exercice correspond typiquement au travail qu’on pourrait donner à des étudiants de licence 2 ou 3 ayant suivi des cours de statistique, cartographie et SIG.

On propose de réaliser une analyse de la distribution de l’Indicateur de Développement humain des régions de Côte d’Ivoire entre 2002 et 2022 à l’aide des données accessibles sur le site Global Data Lab . Le programme est rédigé de telle sorte qu’on puisse facilement refaire l’exercice en prenant d’autres dates (ex. Evolution entre 2012 et 2022) ou d’autres pays d’Afrique de l’Ouest.

Code

library(sf)

	gdlcode	Region	pop2002	pop2022	hdi2002	hdi2022
20	CIVr101	Abidjan	3757	5984	0.581	0.569
21	CIVr102	Yamoussoukro	154	245	0.636	0.619
22	CIVr103	Bas Sassandra	729	1160	0.429	0.520
23	CIVr104	Comoe	1031	1642	0.415	0.582
24	CIVr105	Denguele	525	836	0.353	0.440
25	CIVr106	Goh-Djiboua	944	1504	0.448	0.551
26	CIVr107	Lacs	1513	2409	0.418	0.542
27	CIVr108	Lagunes	1410	2245	0.479	0.566
28	CIVr109	Montagnes	1880	2994	0.400	0.492
29	CIVr110	Sassandra-Marahoue	2224	3541	0.427	0.512
30	CIVr111	Savanes	864	1376	0.423	0.471
31	CIVr112	Vallee du Bandama	716	1140	0.465	0.554
32	CIVr113	Woroba	791	1259	0.243	0.423
33	CIVr114	Zanzan	1149	1829	0.348	0.491

(2) STATISTIQUE UNIVARIEE

Paramètres principaux

• Consigne : Etudiez l’évolution de l’IDH entre 2002 et 2022 en vous servant de paramètres principaux (valeurs centrales, paramètres de dispersion). Puis établissez deux histogrammes permettant de visualiser l’évolution.

Code

# sélectionne les variables
sel <- don[,c("hdi2002","hdi2022")]

# Tableau standard
quant<-apply(sel,2,quantile)
moy<-apply(sel,2,mean)
ect<-apply(sel,2,sd)
cv<-100*ect/moy
tab<-rbind(quant,moy,ect,cv)
row.names(tab) <-c("Minimum","Q1","Médiane","Q3","Maximum","Moyenne","Ecart-type", "C.V. (%)")

kable(tab, caption="Paramètres principaux", digits =2,
      col.names = c("Situation en 2002","Situation en 2022"),
      )

Paramètres principaux

	Situation en 2002	Situation en 2022
Minimum	0.24	0.42
Q1	0.40	0.49
Médiane	0.42	0.53
Q3	0.46	0.56
Maximum	0.64	0.62
Moyenne	0.43	0.52
Ecart-type	0.10	0.06
C.V. (%)	21.99	10.65

Histogrammes

• Consigne : Etablissez deux histogrammes permettant de visualiser la forme de la distribution de l’IDH en 2002 et en 2022

Code

par(mfrow=c(2,1))

mintot <-min(c(sel$hdi2002, sel$hdi2022))
maxtot <-max(c(sel$hdi2002, sel$hdi2022))

# Histogramme
hist(sel$hdi2002,
     breaks=quantile(sel$hdi2002),
     xlim=c(mintot,maxtot),
     col="gray80",
     main= "Situation en 2002",
     xlab = "IDH régional",
     ylab = "Fréquence moyenne")
rug(sel$hdi2002, col="black", lwd=2)
lines(density(sel$hdi2002), lty=3,lwd=2)

hist(sel$hdi2022,
     breaks=quantile(sel$hdi2022),
     xlim=c(mintot,maxtot),
     col="gray80",
     main= "Situation en 2022",
     xlab = "IDH régional",
     ylab = "Fréquence moyenne")
rug(sel$hdi2022, col="black", lwd=2)
lines(density(sel$hdi2022), lty=3,lwd=2)

(3) DONNEES GEOMETRIQUES

Acquision

• Consigne : Après avoir chargé le fonds de carte affichez le contour des régions et le code des unités.

Code

# ne conserve que le code et le nom
map <- map[,c("gdlcode","Region","geometry")]


# Affichage du fonds de carte
par(mar=c(0,0,3,0))
plot(map$geometry, 
     col="gray90",
     main = "Code des unités spatiales de la zone d'étude")

# Ajout du code des unités spatiales
coo<-st_coordinates(st_centroid(map))
text(coo, map$gdlcode, cex=0.5,col="black",)

(4) CARTOGRAPHIE THEMATIQUE

Cartes de stock

• Consigne : Réalisez deux cartes de stock décrivant le nombre d’habitants en 2002 et 2022 Vous utiliserez la même échelle de taille pour rendre les deux cartes comparables.

Code

library(mapsf)
map<-map[,c("gdlcode","geometry")]
map_don <- merge(map, don, by="gdlcode")

maxequ<-max(don$pop012,don$pop2022)

par(mfrow=c(1,2))
mf_map(map_don$geometry, col="white")
mf_map(map_don, type="prop", var="pop2002",
       val_max = maxequ, inches=0.1, col="gray20", 
       leg_title = "Population",)
mf_layout(title="2012",frame = T, credits = "Source : INS Tunisie")

mf_map(map_don$geometry, col="white")
mf_map(map_don, type="prop", var="pop2022",
       val_max = maxequ, inches=0.1, col="gray20",
       leg_title = "Population")
mf_layout(title="2022",frame = T, credits = "Source : INS Tunisie")

Cartes de ratio (choroplèthes)

• Consigne : Réalisez deux cartes de taux décrivant le niveau de l’IDH en 2002 et 2022. Pour les rendre comparables vous utiliserez dans chaque carte une partition en quintiles (5 classes d’effectifs égaux)

Code

library(mapsf)
map_don <- merge(map, don, by="gdlcode")
maxequ<-max(don$hdi2002,don$hdi2022)

par(mfrow=c(1,2))
mf_map(map_don, type="choro", var="hdi2002",
       breaks = "quantile",nbreaks = 5, pal ="Grays",
       leg_title = "IDH",leg_val_rnd = 2)
mf_layout(title="2002",frame = T, credits = "Source : Global Data Lab")

mf_map(map_don, type="choro", var="hdi2022",
       breaks = "quantile",nbreaks = 5, pal ="Grays",
       leg_title = "IDH",leg_val_rnd = 2)
mf_layout(title="2022",frame = T, credits = "Source : Global Data Lab")

(5) STATISTIQUES BIVARIEES

Nuage de points

• Consigne : Tracez un nuage de point montrant l’évolution de l’indicateur entre les deux dates.

Code

# prépration de l'analyse
gdlcode<-don$gdlcode
nom<-don$Region
X<-don$hdi2002
Y<-don$hdi2022
tab<-data.frame(gdlcode,nom,X,Y)

# Diagramme
plot(tab$X,tab$Y, 
     pch=20,
     cex=0.8,
     col="red",
     main = "Evolution de l'IDH",
     xlab="IDH 2002",
     ylab ="IDH 2022")
text(tab$X,tab$Y,tab$nom, 
     pos=2,
     cex=0.5,
     col="blue")

Analyse de la corrélation

• Consigne : calculez les coefficients de corrélation de Pearson et Spearman et testez leur sgnificativité.

Code

cor.test(X,Y, method="pearson")

    Pearson's product-moment correlation

data:  X and Y
t = 5.4499, df = 12, p-value = 0.0001477
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.5674649 0.9494016
sample estimates:
     cor 
0.843945 

Code

cor.test(X,Y, method="spearman")

    Spearman's rank correlation rho

data:  X and Y
S = 102, p-value = 0.001742
alternative hypothesis: true rho is not equal to 0
sample estimates:
      rho 
0.7758242 

Droite de régression

• Consigne : calculez l’equation de la droite de régression et tracez- là sur le graphique.

Code

modreg <- lm(Y~X)
summary(modreg)

Call:
lm(formula = Y ~ X)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.047666 -0.013634 -0.003286  0.018386  0.067288 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.30960    0.04016   7.709 5.48e-06 ***
X            0.49425    0.09069   5.450 0.000148 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.03115 on 12 degrees of freedom
Multiple R-squared:  0.7122,    Adjusted R-squared:  0.6883 
F-statistic:  29.7 on 1 and 12 DF,  p-value: 0.0001477

Code

# Diagramme
plot(tab$X,tab$Y, 
     pch=20,
     cex=0.8,
     col="red",
     main = "Evolution de l'IDH",
     xlab="IDH 2002",
     ylab ="IDH 2022")
text(tab$X,tab$Y,tab$nom, 
     pos=2,
     cex=0.5,
     col="blue")

abline(modreg,col="black",lwd=1)

Analyse des résidus

• Consigne : Calculez les valeurs théoriques prévus par le modèle de régression et les résidus. Affichez le tableau correspondant après l’avoir trié par ordre de résidus croissants.

Code

tab$Y_est <- modreg$fitted.values
tab$Y_res <- modreg$residuals
tab<-tab[order(tab$Y_res),]
kable(tab, digits=3)

	gdlcode	nom	X	Y	Y_est	Y_res
11	CIVr111	Savanes	0.423	0.471	0.519	-0.048
5	CIVr105	Denguele	0.353	0.440	0.484	-0.044
1	CIVr101	Abidjan	0.581	0.569	0.597	-0.028
9	CIVr109	Montagnes	0.400	0.492	0.507	-0.015
10	CIVr110	Sassandra-Marahoue	0.427	0.512	0.521	-0.009
13	CIVr113	Woroba	0.243	0.423	0.430	-0.007
2	CIVr102	Yamoussoukro	0.636	0.619	0.624	-0.005
3	CIVr103	Bas Sassandra	0.429	0.520	0.522	-0.002
14	CIVr114	Zanzan	0.348	0.491	0.482	0.009
12	CIVr112	Vallee du Bandama	0.465	0.554	0.539	0.015
8	CIVr108	Lagunes	0.479	0.566	0.546	0.020
6	CIVr106	Goh-Djiboua	0.448	0.551	0.531	0.020
7	CIVr107	Lacs	0.418	0.542	0.516	0.026
4	CIVr104	Comoe	0.415	0.582	0.515	0.067

Cartographie des résidus

• Consigne : Cartographiez les résidus après les avoir standardisés.

Code

library(mapsf)

# Standardisation des résidus
tab$Y_res_std<-tab$Y_res/sd(tab$Y_res)

# Jointure avec la carte
map<-map[,c("gdlcode","geometry")]
map_reg <- merge(map, tab, by="gdlcode")

# Choix de la palette et des classes
library(RColorBrewer)
mypal<-brewer.pal(n = 6, name = "RdYlBu")
mybreaks = c(-10, -2,-1,0,1,2,10)

mf_map(map_reg, type="choro", var="Y_res_std",
       pal = mypal, breaks=mybreaks,
       leg_title = "Résidus standardisés",leg_val_rnd = 1)
mf_layout(title="Ecarts à la tendance 2002-2022",frame = T, credits = "Source : Global Data Lab")

(6) ANALYSE SPATIALE

Distance au chef_lieu

• Consigne : ajoutez au tableau de données une colonne correspondant à la distance en km à la ville principale du pays (ici Abidjan) et faites en une cartographie en prenant comme bornes de classes 0, 100, 200, 300, 400, 500, 1000 km.

Code

map_don$cap<-0
map_don$cap[1]<-1
cap<-map_don[map_don$cap==1,]
map_don$dist<-1+as.numeric(st_distance(st_centroid(map_don),st_centroid(cap)))/1000

# Choix de la palette et des classes
mypal<-brewer.pal(n = 7, name = "Greys")
mybreaks = c(0,100,200,300,400,500, 1000)

mf_map(map_don, type="choro", var="dist",
       pal = mypal, breaks=mybreaks,
       leg_title = "en km",leg_val_rnd = 0)
mf_layout(title="Distance à vol d'oiseau au chef-lieu",frame = T, credits = "Source : INS Tunisie")

Relation entre équipement et distance au chef-lieu

• Consigne : Déterminez le modèle qui décrit le mieux la relation entre l’IDH (Y) et la distance à lamétropole économique (X) en 2022

Code

X<-map_don$dist
Y<-map_don$hdi2022

par(mfrow=c(2,2))
plot(X,Y,main="Modèle arithmétique Y = a.X+b", sub = round(cor(X,Y),3))
plot(Y, log(X), main = "Modèle logarithmique Y = a.log(X)+b",sub = round(cor(log(X),Y),3))
plot(log(Y),X, main = "Modèle exponentiel log(Y) = a.X+b",sub = round(cor(X,log(Y)),3))
plot(log(X),log(Y), main = "Modèle puissance log(Y) = a.log(X)+b",sub = round(cor(log(X),log(Y)),3))